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Abstract 

We  model  electromagnetic  interrogation  of  a  polyurethane  foam  using  the  TE  mode  of  the  2D 
Maxwell’s  equations  reduced  to  the  wave  equation  for  a  fixed  frequency  in  the  THz  regime.  The  foam 
block  target  contains  knit  lines  which  are  modeled  by  modifying  the  speed  of  propagation,  i.e.,  by  altering 
the  index  of  refraction.  We  describe  our  efforts  to  estimate  the  dielectric  constant  in  the  knit  lines,  as 
well  as  in  the  surrounding  foam,  by  use  of  the  classical  Clausius-Mossotti  equation,  assuming  only  a 
change  in  density.  We  compare  the  numerical  simulations  accounting  for  knit  lines  to  those  in  which  knit 
lines  are  neglected,  each  in  the  context  of  modeling  reflections  of  plane  waves  in  foam  with  voids. 


1  Introduction 

The  problem  we  consider  is  the  scattering  of  a  THz  plane  wave  in  the  possible  presence  of  thin  knit  lines  (i.e., 
layers  of  increased  density)  and  voids  (i.e.,  pockets  of  no  density)  inside  a  block  of  low  density  polyurethane 
foam,  similar  to  BX-250,  which  was  used  on  the  Space  Shuttle  Columbia.  The  detection  of  voids  inside 
the  Sprayed  on  Foam  Insulation  (SOFI)  belonging  to  the  Thermal  Protection  System  (TPS)  is  of  critical 
importance  to  the  NASA  Return  to  Flight  effort.  In  initial  efforts,  THz  frequency  waves  have  been  shown 
to  be  particularly  useful  in  foam  interrogation  [8].  However,  the  modeling  of  and  data  interpretation  for  the 
propagation  of  a  THz  pulse  inside  of  a  material  which  exhibits  heterogeneous  microstructures  of  sizes  that 
are  on  the  order  of  the  wave  length  of  the  interrogating  field  is  not  straight-forward.  In  addition,  there  is 
presently  a  paucity  of  data  on  the  dielectric  properties  of  low  density  foam  in  the  THz  regime.  The  work  in 
[9]  begins  to  remedy  this  deficiency. 

Previous  efforts  in  THz  interrogation  of  SOFI  generally  involved  limited  signal  processing  techniques 
such  as  peak-to-peak  intensity  ratio  detection  [3]  and  time  of  flight  methods  to  determine  the  existence 
of  material  variations  (see  for  example  [10]).  Such  approaches  do  not  take  advantage  of  very  much  of  the 
information  potentially  contained  in  the  reflected  signal.  A  physics-based  model  could  be  used  to  more 
accurately  ascertain  the  geometric  properties  of  an  anomaly,  such  as  size  and  depth,  as  well  as  to  distinguish 
between  a  variation  due  to  the  presence  of  a  void  and  a  normal  variation  due  the  presence  of  knit  lines.  The 
corresponding  information-rich  experimental  data  are  difficult  to  obtain  since  the  amplitudes  of  reflections 
from  low  density  materials  are  very  low  to  immeasurable  using  currently  available  power  sources  at  the  THz 
frequency  (e.g.,  see  Figure  2  of  [7]).  The  theory  and  computations  herein,  therefore,  are  in  the  nature  a 
“proof  of  concept”  and  will  hopefully  serve  as  added  justification  and  motivation  for  the  development  of 
more  powerful  generation  devices. 

Models  utilizing  polarization  mechanisms  have  been  previously  investigated  [1,  2],  but  matching  the 
simulations  to  actual  composite  material  experimental  data  has  not  yet  been  completely  successful.  Moreover, 
these  models  were  expressed  in  only  one  spatial  dimension  and  thus  they  do  not  allow  for  non-normally 
incident  angles  or  non-coplanar  interfaces.  Previous  models  also  did  not  explicitly  account  for  the  effect  of 
knit  lines,  using  instead  an  effective  dielectric  constant  computed  from  observed  time  of  flight  measurements. 
The  current  effort  employs  a  two  dimensional  model  of  the  propagation  of  a  THz  pulse  through  a  medium 
with  curved  knit  lines  and  arbitrarily  shaped  voids.  Because  an  appropriate  model  for  the  dispersion  in  this 
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type  of  media  has  not  yet  been  determined,  we  here  neglect  this  phenomena  in  favor  of  focusing  our  attention 
on  the  reflections  and  refractions  at  the  interfaces.  We  note  that  while  a  dispersion  mechanism  may  later  be 
coupled  to  this  model,  the  single  dominant  frequency  aspect  of  this  problem  will  likely  lead  to  similar  results 
in  both  modeling  approaches.  In  particular,  insomuch  as  dispersion  models  treat  the  index  of  refraction  as 
frequency  dependent,  we  may  assume  that  the  constant  index  of  refraction  used  herein  is  fixed  at  the  value 
associated  with  the  dominant  frequency  mode  of  the  interrogating  pulse. 

2  Model 

For  our  domain  we  choose  a  square  region  (0  <  x  <  b,  0  <  y  <  b)  consisting  of  a  low  density  material  with 
possible  layers  of  higher  density  and  pockets  which  are  modelled  as  having  zero  density, i.e.,  voids.  Figure  1 
depicts  a  simulation  with  a  schematic  of  two  knit  lines  (represented  by  dashed  lines)  1  mm  from  each  other, 
each  parallel  to  an  approaching  plane  wave,  and  perpendicular  to  the  direction  of  propagation.  The  far  right 
boundary  (y  =  b )  is  assumed  to  be  metallic  and  therefore  supra-conducting,  thus  simulating  the  aluminum 
backing  of  the  SOFI  on  the  shuttle  external  tank. 
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Figure  1:  Simulation  of  plane  wave  approaching  parallel  knit  lines. 

We  combine  the  TE  mode  of  the  two  dimensional  Maxwell’s  equations  into  one  equation 

e(x)^(M-)+V^VE(M-))  =~f  (1) 

where  x  =  (x,y),  and  e(x)  and  fi(x)  are  the  spatially  dependent  dielectric  permittivity  and  magnetic  perme- 
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ability,  respectively.  The  corresponding  speed  of  propagation  is 


° ^  n(x)  |  e(f)/i(f)’ 

where  Co  is  the  speed  in  a  vacuum  and  n(x)  is  the  index  of  refraction. 

For  our  source  current  Js  we  wish  to  simulate  a  windowed  pulse,  in  this  case  a  pulse  that  is  allowed 
to  oscillate  for  one  half  of  one  period  and  then  is  truncated.  Although  generators  produce  a  curved,  often 
spherical  wave,  we  assume  the  target  is  sufficiently  far  from  the  generator  (approximately  6  in)  so  that  the 
wave  is  essentially  planar  when  it  reaches  our  domain  of  interest.  We  originate  the  pulse  at  x  =  0,  the  left 
edge  of  our  computational  domain,  and  we  model  it  in  space  as  a  delta  distribution  centered  at  x  =  0.  In 
order  to  have  a  smooth  (in  time)  source  we  use  a  function  of  the  form 

J,{t,x)  =  <Kz)e“((t“to)/to}\  (2) 

where  to  =  £// 4  when  tf  is  the  period  of  the  interrogating  pulse.  For  example,  if  the  frequency  is  /  =  .2THz, 
then  t/  =  l//  =  .5x  10-11s.  A  reasonable  value  for  the  exponent  is  7  =  4. 


2.1  Boundary /Initial  Conditions 

Our  domain  is  defined  to  be  the  region  x  =  (x,y)  G  [0,  b]  x  [0,  b\.  Thus  to  model  a  metallic  backing  behind 
the  foam  at  x  =  6,  we  use  reflecting  (Dirichlet)  boundary  conditions 

[E\x=b  =  0. 


In  order  to  have  a  finite  computational  domain,  we  impose  first  order  absorbing  boundary  conditions  at 
x  =  0;  these  are  modeled  as 


'dE 

dt 


=  0. 

x=0 


With  these  boundary  conditions,  ideally  a  normally  incident  signal  passes  out  of  the  computational  domain, 
and  does  not  return,  i.e.,  we  force  it  to  be  absorbed  by  the  boundary.  Note  that  for  signals  that  are  incident 
at  an  angle,  some  reflection  occurs.  Lastly,  to  allow  for  propagation  along  the  top  and  bottom  boundaries 
(y  =  0  and  y  =  b),  we  use  insulating  boundary  conditions 


We  assume  zero  initial  conditions  so  that 


dE 


l  °y  \  y= 0 

dE 


[dy\v= 


=  0 

=  0. 


y=b 


E(0,X)  =  0 
E(  0,£)  =  0. 


2.2  Modeling  Knit  Lines 

To  model  the  speed  of  wave  propagation  in  the  knit  lines  versus  that  in  the  material  surrounding  them,  we 
must  distinguish  between  the  respective  indices  of  refraction.  However,  we  can  currently  only  measure  the 
effective  index  of  refraction  of  the  composite  material,  for  example,  by  computing  the  “time  of  flight”  in 
experiments.  Thus  we  need  to  relate  these  three  indices  to  each  other  in  order  to  have  accurate  estimates 
for  the  propagation  speed  for  use  in  simulations. 

The  effective  index  of  refraction  ne  can  be  estimated  via  the  Clausius-Mossotti  equation  (see  [6]),  by 
assuming  the  total  polarizability  is  the  weighted  sum  of  the  two  polarizabilities  in  each  part  of  the  material. 
In  particular 

n2e~l  =  P  n\- l  p  n\-  1 
n\  +  2  pi  n\  +  2  p2  +  2  ’ 
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where  pi  and  rti  are  the  densities  and  index  of  refraction,  respectively,  in  each  part  of  the  material,  where 
part  1  corresponds  to  the  low  density  region,  and  part  2  corresponds  to  the  knit  lines.  The  value  p  above  is 
the  total  density  given  by 

p  m  vp2  +  (1  -  v)p\  (4) 

if  we  assume  that  the  knit  lines  comprise  a  certain  volume  fraction  v  of  the  foam.  Let  the  knit  line  density 
P2  be  some  constant  multiple  of  pi  representing  an  increased  density  in  the  knit  lines,  i.e., 


p2  —  P  Pi  • 


(5) 


Substituting  (4)  and  (5)  into  (3)  we  obtain 


nl+2 


(WJ+1 


0 


n\  —  1 
+  2 


rcj-M 

P{n2  +  2)y 


(6) 


If  we  further  assume  that  the  polarizability  of  the  knit  lines  is  equal  to  that  of  the  low  density  region  (this 
is  in  recognition  that  polarizability  is  a  molecular  characteristic  and  the  molecules  in  the  high  density  knit 
lines  are  not  assumed  to  be  significantly  different  than  those  in  the  low  density  regions;  only  the  number  of 
molecules  is  different),  then  we  also  have  the  following 


n\  —  1  nl  —  1 
n\  +  2  =  P(n%  +  2)  ’ 


and  therefore, 


^  =  2(„/3  +  l-„4^. 


n ; 


(8) 


Thus,  if  ne  is  estimated  via  experiments,  n\  can  be  determined  using  equation  (8)  with  reasonable  values  of 
v  and  yd,  and  n 2  can  in  turn  be  calculated  with  equation  (7). 

Experiments  have  suggested,  via  time-of-flight  measurements,  a  value  for  the  effective  index  of  refraction 
of  ne  =  1.03225  ±  0.001.  The  volume  fraction  v  can  be  estimated  by  noting  the  thickness  of  each  knit  line 
divided  by  the  period  in  which  the  knit  lines  occur.  For  example,  0.5  mm  knit  lines  in  each  0.5  cm  of  foam 
corresponds  to  v  =  .1.  The  compression  factor  (3  is  more  difficult  to  measure,  but  experiments  can  be  done 
to  weigh  samples  with  varying  concentrations  of  knit  lines  to  estimate  the  increase  in  density.  An  initial 
estimate  based  on  pictures  of  SOFI  under  20X  magnification  [4]  is  j3  =  2.5.  Using  the  values  yd  =  2.5  and 
v  =  .1,  we  can  estimate  the  index  of  refraction  in  each  part  of  the  foam  to  be 


ni  =  1.01398  (9) 

n2  =  1.03507.  (10) 


If  we  compute  the  decrease  in  speed  of  a  THz  pulse  due  to  the  presence  of  knit  lines  by 

1  —  Til 

P  =  - > 

ne 

then  (9)  would  correspond  to  an  observed  decrease  in  speed  of  1.77%.  Laboratory  experiments  have  suggested 
that  the  decrease  in  speed  is  around  2%. 


3  Results 

We  briefly  describe  here  our  techniques  for  solving  the  system  described  in  (1)  with  the  boundary  and  initial 
conditions  outlined  in  Section  2.1.  Simulations  for  the  cases  in  which  knit  lines  are  ignored  are  compared  to 
those  where  knit  lines  are  described  using  the  values  estimated  above. 
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3.1  Numerical  Solution 


We  apply  a  Finite  Element  method  using  standard  linear  two  dimensional  (Q 1)  basis  elements  to  discretize 
the  model  described  by  (1)  in  space.  This  results  in  nine-banded  mass  and  stiffness  matrices,  M  and  S', 
respectively.  We  also  have  a  contribution  from  absorbing  boundaries  which  we  denote  by  B.  Thus  our 
semi-discrete  system  for  the  vector  of  electric  field  values  e  is 

Me  +  Be  +  Se  =  f. 

(Note  that  we  have  absorbed  coefficients  ^  and  1  into  the  definitions  of  M  and  H,  respectively.) 

For  the  temporal  derivatives  we  use  second  order  discretizations  (centered  differences)  for  both  the  first 
and  second  derivatives.  After  collecting  all  terms  involving  the  updated  time  step  into  the  left  side  of  the 
equation,  we  have  the  following  linear  system 


Aen+ i  =  d,  (11) 

where  A  contains  multiples  of  M  and  H,  and  d  depends  on  en  and  en_i,  as  well  as  S  and  /. 

We  ran  numerical  comparisons  using  various  linear  solvers  including  preconditioned  conjugate-gradient 
and  sparse  LU  factorization.  The  fact  that  the  matrix  A  is  stationary  in  time  contributes  to  the  fact  that 
LU  factorization  performed  better  than  any  iterative  method.  However,  the  size  of  the  problem  that  could 
be  addressed  was  severely  limited  by  the  memory  constraints  when  the  LU  factors  exhibit  fill-in.  The 
iterative  method  on  the  other  hand  can  be  formulated  using  a  matrix-free  approach,  thus  freeing  memory 
for  representing  a  larger  solution.  Unfortunately,  the  computation  time  increases  to  an  unacceptable  level. 

Therefore  we  prefer  to  use  a  mass-lumping  approach  where  quadrature  rules  are  applied  to  the  basis 
functions  to  form  mass  and  stiffness  matrices  which  are  diagonal.  This  results  in  an  explicit  linear  system  for 
(11).  This  system  is  simple  to  solve  at  each  time  step  as  it  requires  only  division  by  the  diagonal  elements. 
There  is  an  obvious  loss  in  accuracy  due  to  the  approximate  integration  inherent  in  mass-lumping,  but 
the  increased  efficiency  allows  for  a  finer  discretization  which  can  sometimes  compensate.  In  fact,  in  our 
testing  of  small  sample  problems  there  is  not  a  noticeable  difference  in  accuracy  between  the  three  methods 
(when  the  computational  times  are  comparable).  However,  a  distinct  difference  is  that  for  the  LU  method 
numerical  error  presented  itself  as  oscillations  preceding  the  signal,  whereas  for  the  mass-lumped  problem, 
the  oscillations  followed  the  signal.  Because  the  beginning  of  the  reflection  is  important  to  determining  the 
location  and  composition  of  a  defect,  we  prefer  for  numerical  error  to  trail  the  propagating  wave  and  thus 
chose  the  lumped  mass  system  for  our  simulations. 

3.2  Simulations 

We  perform  numerical  simulations  of  a  plane  wave  propagating  through  a  material  described  by  its  index  of 
refraction  which  determines  the  speed  of  propagation.  We  consider  the  presence  of  a  void  similar  to  what 
is  seen  in  SOFI  when  a  layer  does  not  completely  fill  a  recess  formed  in  the  previous  layer,  thus  causing  a 
pocket  of  air  to  be  trapped.  The  void  is  modeled  by  taking  its  index  of  refraction  to  be  that  of  free  space 
(n0  =  1). 

Figure  2  displays  snapshots  in  time  of  the  propagation  of  a  plane  wave  (the  white  band)  incident  on  a 
void  in  the  material.  Here  we  neglect  the  specific  properties  of  the  knit  lines  by  modeling  the  entire  foam 
block  using  only  the  effective  (“observed”)  index  of  refraction  ne.  The  reflection  from  the  void  is  clearly 
seen  in  the  second  frame.  This  reflection  expands  out  to  form  an  oblong  elliptical  wave  which  eventually 
propagates  back  to  the  antenna  where  the  signal  is  recorded  with  a  receiver. 

The  simulations  of  the  modeling  approach  proposed  in  this  paper,  namely  the  scenario  where  the  knit 
lines  are  specifically  modeled  with  their  own  index  of  refraction  ri2  and  the  surrounding  low  density  regions 
are  described  by  n i,  are  displayed  in  Figure  3.  Reflections  from  the  knit  lines  are  apparent  in  the  first  frame. 
As  in  the  homogeneous  case  above,  the  reflection  from  the  void  is  again  visible  in  the  second  frame.  The 
interacting  reflections  are  clearly  more  complicated  in  this  scenario,  as  is  the  data  collected  by  the  receiver. 
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Figure  2:  Surface  plots  of  solutions  for  the  case  where  the  effective  index  of  refraction  is  used  in  the  low 
density  part  of  the  foam  and  in  the  knit  lines  (i.e.,  the  presence  of  knit  lines  is  ignored). 
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Figures  4  and  5  depict  the  data  collected  over  time  at  the  receiver  for  the  two  different  modelling 
approaches,  respectively.  In  the  main  plot  of  each  figure  the  magnitude  of  the  reflection  relative  to  that  of 
the  interrogating  signal  is  apparent.  The  inset  plot  of  each  figure  presents  a  magnification  of  the  reflection 
from  the  void.  There  is  a  distinct  difference  in  the  structure  of  the  two  reflections.  In  particular,  the  reflection 
from  the  knit  line  is  clearly  visible  in  front  of  the  reflection  from  the  void  in  the  second  case.  Note  that  the 
amplitude  of  the  reflections  in  each  case  is  roughly  equivalent,  while  in  the  latter  case  the  disruption  to  the 
electric  field  occurs  for  a  longer  duration. 


Signal  received  at  center  line  (x=0) 


Figure  4:  Signal  received  at  x  =  0  on  the  center  line  for  the  homogeneous  case  (the  inset  is  a  magnification 
of  the  reflection). 


Signal  received  at  center  line  (x=0) 


Figure  5:  :  Signal  received  at  x  =  0  on  the  center  line  for  the  case  where  the  knit  lines  and  surrounding 
regions  are  modeled  separately  (the  inset  is  a  magnification  of  the  reflection). 


4  Conclusions 


We  have  developed  a  framework  which  accounts  for  the  presence  of  knit  lines  in  modeling  the  electromagnetic 
propagation  of  interrogating  pulses  in  SOFI.  We  were  able  to  compute,  using  classical  electromagnetics 
formulations,  estimates  for  the  index  of  refraction  in  both  the  knit  lines  and  the  surrounding  low-density 
foam  based  on  measured  values  from  time-of-flight  experiments  and  observable  material  properties.  The 
values  were  shown  to  be  consistent  with  the  laboratory  based  estimates  that  a  2%  decrease  in  speed  should 
be  evident  in  the  absence  of  knit  lines. 

In  the  effort  to  detect  flaws  in  low-density  materials  such  as  foam,  highly  accurate  models  are  required  to 
give  simulations  the  precision  necessary  to  distinguish  small  amplitude  reflections  from  noise,  including  that 
from  model  error.  In  particular,  physics-based  models  have  a  foundation  in  theory,  lending  credibility  and 
more  importantly,  reliability.  To  validate  such  models,  it  is  necessary  to  compare  predictions  to  experimental 
data.  This  validation  process  sometimes  highlights  short-comings,  which  require  new  physics  to  be  developed. 

The  work  herein  provides  an  approach  to  enhance  the  accuracy  of  a  model  by  making  it  more  represen¬ 
tative  of  the  material  in  question,  while  at  the  same  time  not  increasing  substantially  the  computational 
complexity  of  the  system  to  be  solved.  The  results  themselves  suggest  that  knit  lines  should  be  taken  into 
consideration  in  any  precise  modeling  effort.  But  the  approach  of  using  the  Clausius-Mossotti  equation  to 
augment  the  representation  of  the  dielectric  constant  can  also  be  applied  more  generally  to  other  models. 
In  particular,  it  is  already  being  used  to  improve  the  performance  of  high-accuracy  GPS  measurement  by 
accounting  for  the  presence  of  water  vapor  in  the  air  [5].  It  is  entirely  possible  that,  rather  than  requiring  a 
complicated  distribution  of  permittivities  to  account  for  uncertainty  due  to  the  fluctuations  in  an  otherwise 
homogeneous  material,  a  simple  distribution  of  densities  may,  through  the  Clausius-Mossotti  equation,  lead 
to  a  more  accurate  model  of  the  variability  in  the  dielectric  parameters. 

While  the  current  formulation  may  be  used  as  a  forward  solution  in  an  inverse  problem  context,  it  is 
likely  that  the  highest  value  will  he  in  its  ability  to  generate  synthetic  data  with  which  to  test  faster  signal 
processing  approaches  to  damage  detection.  Thus  it  can  be  used  either  to  explore  which  shapes  of  voids  are 
the  hardest  to  detect,  or  to  generate  data  for  scenarios  that  are  difficult  or  expensive  to  manufacture. 

Future  directions  for  this  work  include  consideration  of  full  Maxwell’s  equations  with  coupled  polarization 
and/or  scattering  mechanisms  to  account  for  the  attenuation  observed  in  experiments.  Further,  in  an  actual 
experimental  setup,  the  transmitter  is  some  distance  from  the  receiver,  thus  the  plane  wave  enters  the 
medium  at  a  slight  angle.  Thus  the  delta  distribution  input  used  to  define  the  antenna  should  be  modified 
to  he  along  a  slanted  line. 
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